1 


Use of in situ and airborne multiangle data to assess 


2 MODIS- and Landsat-based estimates of surface albedo 

3 

4 

5 Miguel O. Roman a ’*, Charles K. Gatebe a ' 6 , Yanmin Shuai fl C , 

6 Zhuosen Wang d ’ e , Feng Gao^, Jeff Masek and Crystal B. Schaaf " 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 * Corresponding author. Tel.: +1 301 614 5498; fax: +1 301 614 5269 

21 E-mail address: Miguel.Q.Roman@nasa.gov 

22 


“ NASA Goddard Space Flight Center, Greenbelt, Maryland, USA 
b Universities Space Research Association (USRA), Columbia, Maryland, USA 
c Earth Resources Technology Inc., Laurel, Maryland, USA 
d Center for Remote Sensing, Department of Geography and Environment, Boston University, Boston, MA, USA 
e State Key Laboratory of Remote Sensing Science, Center for Remote Sensing and GIS of Geography College, 
Beijing Key Laboratory for Remote Sensing of Environment and Digital Cities, Beijing Normal University, Beijing, China 

f 

USDA-ARS Hydrology and Remote Sensing Laboratory, Beltsville, Maryland, USA 
s EEOS Department, University of Massachusetts Boston, Boston, Massachusetts, USA 


1 


23 Abstract - 


24 The quantification of uncertainty of global surface albedo data and products is a critical part 

25 of producing complete, physically consistent, and decadal land property data records for studying 

26 ecosystem change. A current challenge in validating satellite retrievals of surface albedo is the 

27 ability to overcome the spatial scaling errors that can contribute on the order of 20% disagree- 

28 ment between satellite and field- measured values. Here, we present the results from an uncertain- 

29 ty analysis of MODerate Resolution Imaging Spectroradiometer (MODIS) and Landsat albedo 

30 retrievals, based on collocated comparisons with tower and airborne multiangular measurements 

31 collected at the Atmospheric Radiation Measurement Program’s (ARM) Cloud and Radiation 

32 Testbed (CART) site during the 2007 Cloud and Land Surface Interaction Campaign (CLAS- 

33 IC’07). Using standard error propagation techniques, airborne measurements obtained by 

34 NASA’s Cloud Absorption Radiometer (CAR) were used to quantify the uncertainties associated 

35 with MODIS and Landsat albedos across a broad range of mixed vegetation and structural types. 

36 Initial focus was on evaluating inter-sensor consistency through assessments of temporal stabili- 

37 ty, as well as examining the overall performance of satellite-derived albedos obtained at all diur- 

38 nal solar zenith angles. In general, the accuracy of the MODIS and Landsat albedos remained 

39 under a 10% margin of error in the SW (0.3 - 5.0 pm) domain. However, results reveal a high 

40 degree of variability in the RMSE (root mean square error) and bias of albedos in both the visible 

41 (0.3 - 0.7 pm) and near-infrared (0.3 - 5.0 pm) broadband channels; where, in some cases, re- 

42 trieval uncertainties were found to be in excess of 20%. For the period of CLASIC’07, the prima- 

43 ry factors that contributed to uncertainties in the satellite-derived albedo values include: (1) the 

44 assumption of temporal stability in the retrieval of 500 m MODIS BRDF values over extended 



45 periods of cloud-contaminated observations; and (2) the assumption of spatial and structural un- 

46 iformity at the Land sat (30 m) pixel scale. 

47 1 . Background 

48 A major goal of international Earth observation efforts is the long term monitoring of terre- 

49 strial essential climate variables and the production of consistent land surface radiation parame- 

50 ters for rigorous modeling studies. With the advent of a new generation of multi-sensor data and 

51 products for Land science applications, recent efforts have explored the “MODISization” of na- 

52 dir-looking satellite sensors to obtain high-resolution (30 m) MODIS-driven land surface para- 

53 meters (Gao et al. 2006; Roy et al. 2008). For example, (Shuai et al. 2011) combined both Land- 

54 sat reflectance (Masek et al. 2006) and high quality 500 m MODIS BRDF (Bidirectional Reflec- 

55 tance Distribution Function) parameters (Lucht et al. 2000; Schaaf et al. 2002; Schaaf et al. 

56 2008) to retrieve 30 m resolution estimates of surface albedo. By capturing seasonal trends at the 

57 characteristic scale of vegetation change (~1 ha), these approaches have the potential to improve 

58 our understanding of the climate consequences of global land cover change and ecosystem dis- 

59 turbance (Barnes and Roy 2008; Masek et al. 2008). . 

60 The quantification of uncertainty of global surface albedo data and products from both 

61 MODIS and Landsat satellites is a critical part of producing complete, physically consistent, 

62 global, and decadal land property data records. The MODIS BRDF/albedo standard product, 

63 available globally since 2000 at resolutions from 0.5 to 5 km, has been validated to Committee 

64 on Earth Observation Satellites (CEOS) Stage 2 (i.e., over a widely distributed set of locations 

65 and time periods via several ground-truth and validation efforts). This validation stage is a pre- 



66 requisite for any data product that is used for monitoring change over time (Morisette et al. 

67 2002). The high-quality primary algorithm for the MODIS standard albedo product (MCD43) 

68 has also been shown to produce consistent global quantities over a variety of land surface types 

69 and snow-covered conditions (Jin et al. 2003a; Jin et al. 2003b; Salomon et al. 2006; Liu et al. 

70 2009; Roman et al. 2009; Roman et al. 2010; Wang et al. 2011). On the other hand, the combined 

71 MODIS/Landsat albedo product (hereby termed ‘Landsat albedo’), which is based on per-class 

72 MODIS BRDF shapes based on uniform land cover characteristics, has been shown to provide a 

73 more detailed landscape texture and achieve good agreement with in-situ data over a limited 

74 number of field stations (Shuai et al. 2011). Additional assessments over a wide range of spatial 

75 (from 10s of meters to 5-30 km) and temporal scales (from daily to monthly) are nonetheless re- 

76 quired to accurately provide end users with a pixel- specific measure of product uncertainty - 

77 both in terms of retrieval quality (e.g. given a limited number of cloud-free satellite observations) 

78 and their ability to capture albedo trends under conditions of seasonal and/or rapid surface 

79 change. 

80 A current challenge in validating satellite albedo retrievals is the ability to overcome the spa- 

81 tial scaling errors that contribute disagreement between satellite and field- measured values, 

82 which can be on the order of 20% (Jin et al. 2003b; Salomon et al. 2006; Liu et al. 2009; Roman 

83 et al. 2010). Recent studies have acquired measurements atop tall (> 400 m) towers to properly 

84 “scale-up” to satellite measurements (Augustine et al. 2005; Roman et al. 2009). Other efforts 

85 have used high resolution imagery to consider the spatial representativeness of the tower obser- 

86 vation footprint to the MODIS pixel (Susaki et al. 2007; Roman et al. 2009). While these me- 

87 thods provide a good means by which direct “point-to-pixel” assessments can be performed with 



88 high confidence; they present their own set challenges (e.g., in the United States, instruments 

89 atop tall towers cannot be left operating year-round, due to heavy icing and bad weather). On ac- 

90 count of the uncertainties arising from direct comparison between sparsely sampled in situ mea- 

91 surements and their corresponding satellite products, a formal assessment has yet to be carried 

92 out to characterize the ability of the MODIS and Landsat data to capture diurnal trends in albedo 

93 across spatially heterogeneous environments. To address these issues, we present the results from 

94 an uncertainty quantification of MODIS and Landsat albedo retrievals based on collocated com- 

95 parisons with tower and airborne measurements. For the airborne datasets, we have employed the 

96 retrieval scheme presented in Roman et al. (2011a), which follows the operational sequence used 

97 to retrieve the MODIS surface reflectance and BRDF/albedo products, based on high-quality 

98 multiangular reflectance measurements obtained by NASA’s Cloud Absorption Radiometer 

99 (CAR) (King et al. 1986; Gatebe et al. 2003). This study focuses on CAR retrievals obtained 

100 over the Atmospheric Radiation Measurement Program’s (ARM) Cloud and Radiation Testbed 

101 (CART) site during the 2007 Cloud and Land Surface Interaction Campaign (CLASIC’07) (Bin- 

102 dlish et al. 2009; Heathman et al. 2009). 

103 2. Albedo retrieval strategy 

104 In this section, we briefly review the albedo retrieval methods used by the MODIS, Landsat, 

105 and CAR instruments, and assess the calibration performance of the CAR spectral channels dur- 

106 ing the period of CLASIC’07. Readers are referred to Sections 2 and 3 in Roman et al. (201 la) 

107 for detailed descriptions of the CLASIC’07 experiment (including retrieval of CAR and MODIS 



108 BRDF/albedo datasets); and Section 2 in Shuai et al. (201 1) for a complete description of the 

109 Landsat albedo retrieval strategy. 

110 2.1 Instantaneous albedos from CAR, MODIS, and Landsat 

111 The CAR, MODIS, and Landsat albedo retrieval schemes employ the BRDF kernel model 

112 parameters from the reciprocal version of the semiempirical RossThick-LiSparse model 

1 13 (RTLSR) (Wanner et al. 1995; 1997; Lucht et al. 2000): 

114 R a (Q v , ) = f isoA + f voLA K vol (Q v , ) + f ge0yA K geo (Q v , Q s ) (1) 

115 Here, Q v and Q v are the viewing and solar geometries, which are each defined by zenith and azi- 

1 16 muthal angles ( K vo i is the coefficient for the RossThick volume scattering kernel (Ross 

1 17 1981); K geo is the coefficient of the LiSparse-Reciprocal geometric scattering kernel (Li and 

118 Strahler 1992); and/ r ,A are the RTLSR kernel weights x in waveband A with limits [ A miii , A rnax J 

119 (Wanner et al. 1995; Lucht et al. 2000). The RTLSR kernel weights are then used to compute in- 

120 trinsic surface albedos (i.e., black sky albedo for direct beam at local solar noon and white sky 

121 albedo for isotropic diffuse radiation) (Martonchik et al. 2000; Schaepman-Strub et al. 2006): 

I 2 n 1 

Ra(Q,.) = - f d(j) v \R A (Cl v ,Q. i )n v dn Y 

122 n- { \ i (2) 

= fisoA + fvol\ K mlA (Q ) + f geo K K « eoA (^ s ) 

123 

2n 1 

Ra = — [ dfa \RA{0. i )n i d 

124 n i 1 (3) 

= fisoA + fvol A % vo1 A + f geo A % S e0 A 



1 25 where, R A (fi v ,£I ) = nBRDF A (unitless), is the ratio of the surface BRDF to that of a perfect 

126 Lambertian reflector, which can be approximated by measurement over some (small) finite angle 

127 with diffuse illumination and multiple interaction effects accounted for or assumed zero (Lyapus- 

128 tin and Privette 1999). Subscripts v and i denote the upper ‘viewing’ and ‘incident’ hemispheres. 

129 Ra(Q s ) is the black-sky albedo, Ra is white-sky albedo, K vol (Q v ) and K {, Q v ) are the direc- 

130 tional-hemispherical integrals, and K vo i( Q v ) and K geo ( Q v ) are the bihemispherical integrals of 

131 K vol and K geo . Other terms in Eq. (2) and (3) are: 

132 /u y = cos(.9 v ); y = vori (4) 

133 To accurately compare these intrinsic quantities against ground-based albedos, the black-sky and 

134 white-sky albedos must be combined as a function of solar geometry and atmospheric state to 

135 compute instantaneous albedo under assumptions of isotropic diffuse illumination: 

fuo A + (1 ' - A) A ) [LlA K*ol (Q s ) + f geoA K geo (Q f ) ] + 

136 A A/ „(Q i )= = = n (5) 

Aa fvolA K vol ( A ) + f geo A K geo (Q s ) 

137 where D 0a (unitless) is the proportion of diffuse illumination for an absorbing lower boundary 

138 (Lewis and Barnsley 1994; Lucht et al. 2000). The MODIS BRDF shape derived from clear-sky 

139 observations can then be used to derive albedo values in all sky conditions (Liu et al. 2009). 

140 Most recently, the computation of MODIS instantaneous albedos was updated to account for the 

141 effects of multiple scattering and anisotropic diffuse illumination (Roman et al. 2010): 



142 A A (n s )*f isoA +f wlA K" volA ) + fgeoA % geoA(Cl s ) 


(6) 


143 K" xa (Q s ) = (1 - D 0 A )K x (Cl s ) + D qa K'xa (7) 

i 
0 

1 

0 

146 where, V^lfl) is the normalized sky radiance distribution under an absorbing lower boundary 

147 and K xA is the weighted bihemispherical integral of K ^ ( 11.. 11 ^ ) ( \\1 1 e re x — vol or g£?o). 

148 Intrinsic albedo quantities derived from RTLSR BRDF model inversions can then be combined 

149 with in-situ estimates of cloud fraction (< 0.6), 550 nm aerosol optical depth (AOD), solar zenith 

150 angle (SZA), and Z) 0A to compute clear-sky instantaneous albedos fromMODIS, Landsat, and 

151 CAR data. 

152 The kernel-driven models employed by the MODIS and Landsat albedo products are also 

153 identified as part of the heritage algorithms used to generate the Visible Infrared Imager Radi- 

154 ometer Suite’s (VIIRS) Land Environmental Data Records (EDRs); which aim to ensure continu- 

155 ity for AVHRR and MODIS observations by providing high temporal resolution and wide area 

156 coverage (Lee et al. 2010). The VIIRS Land EDRs are currently being evaluated by NASA and 

157 NOAA to assess their suitability for operational weather forecasting and long-term climate moni- 

158 taring applications (Roman et al. 2011b). 
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144 K'x a = - J dfa 



159 2.2 Narrowband to Broadband Conversion 

160 Since field-measured albedos are commonly measured as broadband quantities, an equiva- 

161 lent set of broadband albedos for MODIS and Landsat were generated for the UV- Visible (0.3 - 

162 0.7 pm), NIR (0.7 - 5.0 pm), and the entire spectrum of solar radiation ([SW] 0.3 - 5.0 qm), 

163 based on empirical relations between ground-based albedo measurements and satellite observa- 

164 tions - c/., Eqs. (11) and (15) in Liang (2001). Broadband albedos were also derived for CAR 

165 measurements by calculating the ratio of broadband upwelling radiative flux to broadband 

166 downwelling flux (Liang 2001; Liang et al. 2003): 

A max 

jA(fi,)D(n„A)iA 

167 F(n>A*_ = y>,A(n„A,) (10) 

|£>(n,,A)dA 

Amin 

168 Then, CAR narrowband-to-broadband spectral albedo coefficients, c,-, were generated for each 

169 spectral band by determining the downward fluxes (i.e. direct and diffuse) using an library of 30 

170 reflectance spectra of representative land covers in the ARM Southern Great Plains (SGP) region 

171 (Trishenko et al. 2003): 

172 A short = 0. 160a, +0.291a 2 +0. 243a, +0.116a 4 +0.112a 5 +0.081a 7 -0.0015 (11) 

173 A nir =0.039a, +0.504a 2 -0.071a, +0.105a 4 +0.252a 5 +0.069a 6 +0.101a 7 (12) 

174 \ isible = 0.331a, + 0.424a, + 0.246a 4 


( 13 ) 



175 The upward fluxes were directly obtained from the library of 30 SGP reflectance spectra; while 

176 the downward fluxes were obtained by performing multiple MODTRAN®5.1 (Berk et al. 2004) 

177 mns for a broad range of snow-free conditions (i.e., 21 atmospheric visibility values for different 

178 aerosol loadings, 2 atmospheric profiles, and solar zenith angles ranging from 0° - 80° with the 

179 increment of 1°). 

180 2.3 CAR instrument performance during CLASIC’07 

181 During the CLASIC’07 experiment, radiometric calibration of the CAR spectral channels 

182 was made at the NASA Goddard Space Flight Center Radiometric Calibration Facility (GSFC- 

183 RCF) (Butler and Barnes 2003). A description of the calibration scheme, using a series of inte- 

184 grating spheres with diameters of 1.83 m, 1.22 m, and 0.51 m, covering all of the CAR’s spectral 

185 channels, can be found in Gatebe et al. (2007). The conversion from Digital Numbers (DNs) to 

186 Level 1 at-sensor radiances is determined from the instrument’s response for at least two known 

187 radiance levels and then determining the instrument gain (slope) and offset (intercept) for each 

188 wavelength across the sensor band pass. The estimated errors associated with this radiometric 

189 conversion vary from ±1% to ±3% for all spectral channels (Gatebe et al. 2003; Gatebe et al. 

190 2007). Radiometric calibration was performed prior to and after CLASIC’07. In the past, to de- 

191 termine a suitable calibration for a given flight during the experiment, a linear change between 

192 the preflight and postflight calibration was assumed as a function of only the number of flights 

193 flown during an entire campaign. For the CLASIC’07 experiment, however, both the pre- and 

194 post-calibration coefficients were averaged. This was found to be representative of each flight 

195 scenario, and made it easier to account for uncertainties related to calibration, stability, and wa- 



196 velength errors. We note that the calibration ratios, postflight-to-preflight, varied between 0.9691 

197 (at A = 0.472 pm) and 1.1845 (at A = 0.340 pm). 

198 Table 1. Remotely-sensed albedo retrieval scenarios obtained during CLASIC’07 for MODIS, 

199 Landsat, CAR in medium resolution mode (MRM), and CAR in coarse resolution mode (CRM). 


Instrument 

GIFOV (m) 

Scale (m) 

BRDF Retrieval Period 
(DOY) 

Albedo Reconstruction Period 
(DOY) 

MODIS 
Terra/ Aqua 

436 - 1686 

500 

145-193 

153-155; 159-190 

Landsat TM 

30 

30 

153-168 

153-155 

CAR-MRM 

15-45 

30 

175 

153-155 

CAR-CRM 

90 - 360 

250 

175 

159-190 


200 

201 3. Retrieval scenarios during CLASIC’07 

202 Table 1 provides a summary of the individual BRDF retrieval and albedo reconstruction pe- 

203 riods for CAR, MODIS, and Landsat. Note that the retrieval scenarios varied by sensor. For in- 

204 stance, Landsat albedos were reconstructed for a short time period (Day of the year, DOY 153- 

205 155) to better represent the per-class albedo-to-nadir-reflectance (A/N) ratios derived from the 

206 concurrent MODIS acquisition period. Conversely, the CAR albedos are based on two different 

207 modes: a Medium Resolution Mode (CAR-MRM) to match the scale of Landsat data; and a 

208 Coarse Resolution Mode (CAR-CRM) to match the scale of MODIS data (c/.. Fig. 9 in Roman et 

209 al. 2011a). The CAR measurements are based on CLASIC Flight #1928 (DOY 175, 2007). Ac- 

210 cordingly, a longer measurement period was also examined (i.e., DOY 159-190) to evaluate the 

211 ability of the CAR and MODIS data to represent the landscape conditions surrounding the entire 

212 CLASIC’07 period. 


213 




214 


a. Landsat b. CAR-MRM c. CAR-CRM d. MODIS 



White-Sky Albedo (SW: 0.3 - 5.0um) 



0.100 0.133 0.167 0.200 0.233 0.267 0.300 


e. Landsat f. CAR-MRM 



a CAR-CRM 

I 




NDVI based on Broadband White-Sky Albedo 



0.35 0.44 0.52 0.61 0.69 0.77 0.86 

i. IKONOS somrower j 2D Gaussian / 60 "’ Tower 



215 Fig. 1. (a. - h.) 2 km x 2 km subsets illustrating retrievals of white -sky albedo (WSA) and broad- 

216 band NDVI. The Broadband NDVI values are based on the UV-VIS and NIR portions of the solar 

217 spectrum: NDVI = (WSA N ir - WSA V is) / (WSA N ir + WSA V is). (i-) IKONOS true-color scene, re- 

218 trieved on 1 July 2007, denotes the -960 m ground footprint as seen by the downward-facing pyra- 

219 nometer atop the CART site’s 60 m radiation tower, (j.) Spatially-integrated albedos for CAR, 

220 Landsat, and MODIS, were computed based on a 2D Gaussian filter representing the tower foot- 

221 print. 
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Legend: 


■ Landsat BCAR-MRM aCAR-CRM ■ MODIS 


Land Cover Type: 


1) Grass/pasture 

2) Trees/shrubs 

3) Corn/milo 

4) Wheat/stubble 


5) Bare soil/short grass 

6) Bare soil dark 

7) Bare soil medium 

8) Bare soil light 


Fig. 2. Comparisons of intrinsic white-sky albedos derived from Landsat, CAR, and MODIS data 
shown for individual land cover types present across the ARM CART during CLASIC’07. Results 
were partitioned into broadband albedos based on the UV- Visible (0.3 - 0.7 pm), NIR (0.7 - 5.0 pm), 
and the entire spectrum of solar radiation ([SW] 0.3 - 5.0 pm). 


Fig. 1 illustrates several 2 km x 2 km subsets based on CAR, MODIS, and Landsat retrievals 
of SW white-sky albedo (WSA) (Fig. la -d), broadband NDVI (Fig.le -h), as well as a true-color 
IKONOS scene retrieved on 1 July 2007 (Fig.li). The red circle over the IKONOS scene denotes 
the -960 m ground diameter footprint as seen by the downward-facing pyranometer atop the 
CART site’s 60 m radiation tower. Throughout the analysis stage in Section 4, a Gaussian filter 
was applied to the CAR, MODIS, and Landsat data to compute spatially-integrated albedos that 
represented the tower’s ground-projected instantaneous field of view (GIFOV) (Fig. lj). 


235 In order to quantify the temporal consistency of albedo retrievals during the period of 

236 CLASIC’07, estimates of broadband NDVI were generated for each sensor (e.g. Fig. le - h). 

237 NDVI estimates based on broadband white-sky albedos were used to minimize the influence of 

238 variable sun-target- sensor configurations when estimating measurement differences and changes 

239 in vegetation conditions (Huete et al. 2002). Results in Fig. 1 show how areas north of the CART 

240 site and the winter wheat fields just east and west of the site, appear to be more vegetated during 

241 the Landsat overpass period (DOY 154). However, the overall change in NDVI was relatively 

242 small throughout the CLASIC’07 period (DOY 153-190) (i.e., -5.7% for the entire study area 

243 and -8.0% for the area within the CART tower footprint). 

244 The WSA results in Fig. la - Id show some similarities between the CAR-CRM and MOD- 

245 IS WSA data, as well as some differences between the finer resolution Landsat and CAR-MRM 

246 albedos. In particular, the Landsat albedos (Fig. la) could resolve fine-scale spatial features 

247 across the CART site (e.g., small buildings and dirt roads); but were also characterized by higher 

248 ‘within-biome’ variability. This is demonstrated in Fig. 2, which illustrates the overall mean and 

249 standard deviations in white-sky albedo based on individual estimates obtained for each repre- 

250 sentative land cover class identified across the CART site. The ancillary land cover data is based 

25 1 on field surveys, vegetation measurements, and surface characterizations performed during 

252 CLASIC’07 (Roman et al. 2011a). Unlike the MODIS and CAR intrinsic albedos, which are dri- 

253 ven by the anisotropy of each pixel, the Landsat retrieval scheme is based on per-class MODIS 

254 BRDF shapes based on uniform land cover characteristics. This resulted in WSA retrievals of 

255 varying magnitudes, particularly across mixed cover types (e.g., bare soil mixed with short 
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grass). This is in contrast to the relatively lower within-biome variability seen in the CAR and 
MODIS white-sky albedos across all land cover classes and broadband channels. 

Additional quality assurance (QA) checks were performed to assess the consistency of each 
retrieval scheme. For this study, both the Landsat and MODIS intrinsic albedos were based on 
gap-free, quality-enhanced BRDF retrievals that rely on spatial and temporal fitting techniques to 
compensate for missing data and provide an estimate of the surface reflectance anisotropy for 
situations under cloud-contaminated conditions (Zhang 2008). For the period surrounding the 
Landsat date of acquisition, the MODIS retrievals were all based on high-quality “majority” full 
inversion values. Conversely, the data acquisition period surrounding the CAR measurements 
(DOY 175) was impacted by changing weather conditions (e.g., clouds and rainfall events). 

These weather conditions resulted in the majority of MODIS retrievals to be based on lower- 
quality temporally-fitted pixels; particularly, throughout DOY 161-190. Consequently, MODIS 
retrievals provided a close, but not exact, representation of the surface conditions surrounding 
the CAR measurements during the CLASIC’07 experiment. For an in-depth look at the MODIS 
and CAR quality assurance assessment, readers are referred to Section 3 in Roman et al., 

(2011a), where QA summaries based on data from CLASIC Flight #1928 are available. 

Understanding the above mentioned differences in BRDF/albedo retrieval strategies and data 
acquisition periods, it is important to note that the ground-based observations obtained at the 
CART site provide a consistent reference source for all albedo reconstructions (albeit for a small 
area the size of a few MODIS pixels) (c/., Section 4). Furthermore, as we will demonstrate in 



276 Section 5, the CAR data can help reduce the propagation of measurement uncertainty and error 

277 when evaluating the satellite-based retrieval schemes at the individual pixel level. 

278 4. Comparisons to tower-based measurements 

279 We now examine the diurnal performance of instantaneous albedos derived from CAR, 

280 MODIS, and Landsat, based on comparisons against available in situ observations acquired dur- 

281 ing CLASIC’07. Measurements from a downward-facing pyranometer installed on a 60 m radia- 

282 tion tower at the CART site collected albedo and radiation fluxes in the shortwave (SW) (0.3- 2.8 

283 pm) waveband (Fig. li). Two additional instruments, a normal incidence pyrheliometer mounted 

284 on an automatic sun tracker and a shaded pyranometer riding on top of the sun tracker, measured 

285 direct and diffuse solar radiation incident upon the field station ( cf , Fig. 3 in Roman et al., 

286 2011a). Estimates of cloud fraction and aerosol optical depth, as viewed from a skyward-looking 

287 pyranometer and an AERONET sunphotometer (Holben et al. 2001), were also collected. This 

288 measurement scheme follows a strict set of guidelines as established by the International Base- 

289 line Surface Radiation Network (BSRN) (McArthur 2005; WMO 2006). BSRN measurement 

290 protocols are recognized as the international standard for in situ albedo data, with a review 

291 process that includes additional quality assurance (QA) checks (e.g. standard units, naming con- 

292 ventions, and reporting intervals) to maintain consistency within the larger network-wide BSRN 

293 database (Schaaf et al. 2009). 

294 Following the albedo reconstruction periods described in Table 1, results in Figs. 3 and 4 

295 show comparisons between the tower-based albedos and instantaneous albedos derived from 

296 CAR, MODIS, and Landsat data. Measurements of aerosol optical depth (AOD) at 550 nm are 



297 also plotted in Fig. 4 (right-axis); while AOD measurements recorded during DOY 153-155 (Fig. 

298 3) remained low and constant throughout this period (i.e., AOD = 0.1676 ± 0.0590). Results 

299 show the usual “U-shaped” diurnal trend in instantaneous albedo that reaches a minimum value 

300 around local solar noon time. In general, the CAR and MODIS albedos met the absolute accura- 

301 cy requirement of 0.02 units (i.e., within 10% of surface measured values) for instantaneous SW 

302 albedos at SZA < 45° (i.e., between 10.00 and 16.00 local time); with the CAR-CRM albedo data 

303 also performing well at SZA > 45° (i.e., before 10.00 and after 16.00 local time). Conversely, 

304 both MODIS and Landsat consistently underestimated the tower albedos at SZA > 45°. The sta- 

305 tistical results in Tables 2 and 3 show similar negative trends in the biases derived from MODIS 

306 and Landsat albedos. Finally, CAR and MODIS retrievals based on the full expression of instan- 

307 taneous albedo (Eq. 6) showed slight improvements by -0.0065 absolute units over the isotropic 

308 albedo formulation (Eq. 5). As seen for several dates in Fig. 4, the full expression results for 

309 MODIS and CAR are much closer to the daily albedo maxima at SZA = 75°. 
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Fig. 3. Comparisons between instantaneous retrievals of surface albedo (15 -min intervals) derived 
from CAR-MRM (red squares), MODIS (green triangles), Landsat (‘X’ marks), and tower- based 
measurements (blue diamonds) acquired at the CART site throughout a 3-day period surrounding 
the Landsat overpass date (DOY 154, 2007). 

Table 2. Accuracy 3 , (absolute bias) and uncertainty 11 (RMS of absolute error or RMSE) values re- 
sulting from comparisons between between ground-based (CART), airborne (CAR-MRM), and sa- 
tellite-derived (MODIS and Landsat) albedos as illustrated in Fig. 3. The total sample size (n) for 
two solar zenith angle (SZA) ranges is shown. 



10° < SZA < 45° (n = 69) 

45° < SZA < 75° (n = 34) 

DOY 153-155,2007 

CAR-MRM 

MODIS 

Landsat 

CAR-CRM 

MODIS 

Landsat 

Accuracy (Bias) 

0.0044 

-0.0136 

-0.0324 

-0.0266 

-0.0433 

-0.0401 

Uncertainty (RMSE) 

0.0090 

0.0157 

0.0327 

0.0287 

0.0447 

0.0407 


“Accuracy =arithmetic mean (Sensor - To wer) 

Uncertainty: RMS of absolute error = \ arithmctic mean (Sensor - Tower) 2 
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Fig. 4. Comparisons between instantaneous albedos (30-min intervals) derived from CAR-CRM 
and MODIS (using both isotropic and full expressions), and tower-based measurements acquired at 
the CART site throughout a 32-day period surrounding CLASIC Flight #1928 (DOY 175, 2007). 


Table 3. Accuracy and uncertainty values resulting from a 32-day comparison between ground- 
based (CART), airborne (CAR-CRM), and satellite -derived (MODIS) albedos as illustrated in Fig. 


4. Setup is the same as Tal 

lie 2. 

10° < SZA < 45° (n = 289) 

45° < SZA < 75° (n = 193) 

DOY 159-190,2007 

CAR-CRM 

MODIS 

CAR-CRM 

MODIS 

Accuracy (Bias) 

0.0042 

-0.0286 

-0.0096 

-0.0495 

Uncertainty (RMSE) 

0.0082 

0.0296 

0.0184 

0.0526 


5. Regional assessment of MODIS and Landsat albedos 

In the previous section, CAR retrievals were shown to be of sufficient accuracy and consis- 
tency to reproduce the diurnal variations in albedo across the CART site throughout the entire 
period of CLASIC’07. Using the CAR-MRM and CAR-CRM instantaneous albedos as “ground- 
truth”, we now employ standard error propagation techniques (Heuvelink 1998) to quantify the 
uncertainties associated with MODIS and Landsat instantaneous albedos over a mixture of land- 


scapes extending beyond the tower observation footprint at the CART site. 


341 Assuming that the error propagation terms in Eqs. (14 -15) are uncorrelated, the Root-Sum- 

342 of-Squares Error (RSSE) can be used to provide estimates of retrieval uncertainty (absolute 

343 RMSE) and bias, both in an absolute and temporal sense: 



346 where Err(d x J denotes the “in-situ to satellite” error propagation chain based on two component 

347 factors. CAR BRDF/albedo retrievals were then matched to the resolution of MODIS and Land- 

348 sat to minimize errors due to sub-grid scale mismatch and the effects of land surface heterogenei- 

349 ty. Additional checks were also performed to limit the sampling of CAR pixels to the highest 

350 quality “majority” full BRDF inversion values. For CLASIC’07, this resulted in 789 individual 

35 1 samples, each of which was tested following the same albedo reconstruction periods presented in 


352 Section 3. 
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Fig. 5. Distribution of absolute bias (accuracy) and RMSE (uncertainty) for MODIS and Landsat 
instantaneous albedos at UV- Visible, NIR, and SW broadband channels for the CLASIC’07 period 


over the CART site. 


Results were examined by comparing the distribution of biases and RMSEs resulting from 
MODIS and Landsat instantaneous albedos at LSN < SZA < 75°. The histogram plots in Lig. 5 
show a persistent negative bias (-0.03) in SW instantaneous albedos, corresponding to the biases 
recorded for MODIS and Landsat in earlier assessments ( cf , results in Tables 2 and 3). Likewise, 
for the NIR broadband, roughly 19% of Landsat retrievals showed positive biases above the 
standard accuracy limit of +0.072 units. In the VIS broadband, a small fraction of MODIS (23%) 


363 and a large fraction of Landsat retrievals (57%) were also above the standard accuracy limit of 

364 ±0.018 units. To understand the causes of such differences, a collection of ternary diagrams were 

365 created to determine how each sensor performed at the individual pixel level. The diagrams in 

366 Fig. 6 have been arranged such that retrievals located near the top originate from landscapes 

367 dominated by non-photo synthetic (or 'brown') vegetation (i.e., combined land cover (LC) classes 

368 1, 4, and 5 in Fig. 2). Conversely, retrievals located near the bottom- left correspond to areas 

369 dominated by bare soils (i.e., combined LC classes 6, 7, and 8 in Fig. 2), and retrievals located 

370 near the bottom-right correspond to areas dominated by green vegetation (i.e., combined LC 

371 classes 2, and 3 in Fig. 2). Thus, the closer the satellite retrievals are to the center portions of the 

372 ternary diagrams, the more mixed is the landscape. 

373 Results reveal a large degree of variability in the RMSE and bias estimates of MODIS and 

374 Landsat albedos, both between fractional cover types and across broadband channels. In particu- 

375 lar, the MODIS NIR values remained stable across most landscape regimes, with only a few 

376 samples identified above the 20% margin of error for snow-free conditions (Fig. 6e). A synoptic 

377 analysis of the ternary diagrams also suggests that the uncertainties in the NIR broadband are 

378 more likely to propagate into the SW domain (Figs. 6i - 61 ). The MODIS VIS broadband also ap- 

379 peared to capture bare-soil albedo variability (i.e., wet vs. dry areas) with high accuracy (-0.008); 

380 but the biases where moderately larger over mixed landscapes (+0.012) and regions dominated 

381 by non-photo synthetic vegetation (+0.019) (Fig. 6a). Conversely, the Landsat albedos were less 

382 stable in the NIR, with positive biases (+0.05) dominating over mixed landscapes (Fig. 6f). The 

383 same error patterns where seen in the VIS broadband, where more than half of Landsat retrievals 



384 were above the 20% margin of error (Fig. 6d). The latter resulted from varying magnitudes 

385 across mixed cover types and regions dominated by bare soils (Fig. 6b). 
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Fig. 6. Ternary diagrams illustrating the pixel-specific accuracy (absolute bias) and uncertainty 
(absolute RMSE) of MODIS and Landsat instantaneous albedos (LSN < SZA < 75) at UV-Visible (a. 
- d.), NIR (e. - d.), and SW (i. - 1.) broadband channels for the CLASIC’07 period over the CART 
site. For each color scale, green denotes values where the bias or RMSE = 0. For the accuracy dia- 
grams, the lower (blue) and upper (red) limits correspond to retrievals that are at or above a 20% 
margin of error (i.e., relative to in-situ measurements obtained under snow-free conditions). For the 
uncertainty diagrams, the upper (red) limit denotes the same (20%) margin of error. 


395 For the period of CLASIC’07, there were two major factors that contributed to the uncer- 

396 tainties in the satellite-derived albedo values. First, it is unlikely that the assumption of temporal 


397 stability in the retrieval of 500 m MODIS BRDFs could hold together throughout the extended 

398 periods of cloud-contaminated observations experienced during CLASIC’07. Daily records from 

399 land cover surveys performed on DOY 166-173 over landscapes surrounding the CART site con- 

400 firm that several parcels of corn, milo, and winter-wheat were being harvested before full maturi- 

401 ty due to floods experienced along the Salt Fork Arkansas River (located 5 km north of the 

402 CART site). A visual inspection of the 2.4 m IKONOS scene acquired on 1 July (Fig. li) also 

403 con fir ms these events. Thus, it is likely that the negative (-0.03) biases in MODIS SW albedos 

404 are caused by the use of temporally-fitted BRDF shapes that are driven by “majority” full inver- 

405 sion values obtained prior to the early harvesting period (i.e., DOY 145-153, 2007). This may 

406 explain why the uncertainties of MODIS VIS albedos where predominantly above the 20% mar- 

407 gin of error over areas dominated by senescent winter- wheat fields; but remained well under the 

408 10% margin of error over areas dominated by bare soils. 

409 For the Landsat albedos, another major source of uncertainty is the assumption of spatial and 

410 structural uniformity at the Landsat (30 m) pixel scale. In the Landsat albedo retrieval algorithm, 

411 “pure” land cover clusters are identified on a regional basis and then associated with MODIS 

412 anisotropy information through scaling of 500 m BRDF retrievals to 30 m resolution. However, 

413 recent assessments of the CLASIC Flight #1928 dataset have indicated that the use of dominant 

414 archetypal BRDF shapes to describe the anisotropy of heterogeneous pixels may lead to errors on 

415 the order of 0.5% - 6.5% in the retrieved directional reflectance values ( cf , Fig. 10 in Roman et 

416 al., 2011a). This will particularly affect retrievals where heterogeneous conditions are being 

417 lumped into a single land cover class (e.g., bare soil areas not being properly partitioned into dry, 

418 wet, and damp conditions.) As discussed in Shuai et al., (2011) these situations can be addressed 



419 by breaking “pure” land cover clusters into multiple sub-clusters representing different surface 

420 conditions. 

421 6. Conclusions and future recommendations 

422 The diurnal performance of the MODIS and Landsat albedo algorithms (Schaaf et al. 2002; 

423 Shuai et al. 2011) is evaluated using field and airborne measurements coincident with Landsat 

424 TM and multi-date MODIS Terra/ Aqua overpasses. For the broad range of mixed vegetation and 

425 structural types examined during the period of CLASIC’07, the overall accuracy of MODIS and 

426 Landsat SW (0.3 -5.0 pm) albedos is within a 10% margin of error and shows an increasing neg- 

427 ative bias (-0.03) and increased RMSE (0.05) as zenith angle increases compared with the in-situ 

428 measurements. Results also reveal a high degree of variability in the RMSE and bias of MODIS 

429 and Landsat albedos in both the visible (0.3 - 0.7 pm) and near-infrared (0.7 - 5.0 pm) broadband 

430 channels. However, we note that the lack of high-quality “majority” 500m MODIS BRDF pixels 

431 through the experiment hindered the band-dependent quality controls, as outliers were more dif- 

432 ficult to identify. This was particularly the case in the VIS broadband, where cloud contamination 

433 and mixed-pixel contamination are highly likely. Despite such limitations, results obtained indi- 

434 cate that MODIS VIS/NIR albedos are able to capture bare-soil albedo variability (i.e., wet vs. 

435 dry areas) with high accuracy (-0.008). 

436 While recent product development, intercomparison, and validation efforts have focused al- 

437 most entirely on the retrieval of surface albedos for a single SW broadband value, it is important 

438 to note that most numerical prediction models (and global climate and biogeochemical models) 

439 currently in use call for surface energy fluxes and some biophysical variables to be calculated 



440 separately by disentangling broadband albedos into fractional areas of bare soil and vegetation 

441 (Noilhan and Mahfouf 1996; Koster et al. 2000; Ek et al. 2003; Kaptue et al. 2010). It is there- 

442 fore important to continue examining how the accuracies of global albedo products are holding 

443 up in these spectral regimes. Likewise, the uncertainties that may impact satellite-inferred albedo 

444 trends must be assessed and expressed in terms of a reference sensor that can overcome the fore- 

445 told errors due to sub-grid scale mismatch and the effects of land surface heterogeneity. It is thus 

446 critical that continuous, long-term tower measurements of surface albedo and radiation fluxes be 

447 done in concert with intensive airborne measurement campaigns that can focus on addressing 

448 sources of uncertainties at both plot-level (< 90 m) to landscape-level (> 90 m) scales. 

449 It is clear that spatial scale of signal aggregation is very important in the retrieval of mea- 

450 ningful surface radiation properties of vegetated surfaces from multiangle and pseudo-multi- 

45 1 angular remote sensing data. This is forcing experimenters to develop new measurement and va- 

452 lidation protocols for surface BRDF and albedo estimation (Walthall et al. 2000; Schaaf et al. 

453 2009). Ongoing studies combining airborne multiangular measurements from CAR with mea- 

454 surements of terrestrial biomass and ecosystem structure, e.g., NASA’s L-band Digital Beam- 

455 forming Synthetic Aperture Radar (DBSAR) (Rincon et al. 2011) and the Slope Imaging Multi- 

456 polarization Photon-counting Lidar (SIMPL) (Dabney et al. 2010), will provide us with new in- 

457 sights to issues of landscape-level variability and the opportunity to continue examining mixed 

458 pixels from both medium and coarse scale resolution systems. 
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